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We show that the momentum distribution of a nonlinear matter wave suddenly released with a 
finite velocity in a speckle potential converges, after an out-of-equilibrium evolution, to a universal 
Rayleigh-Jeans thermal distribution. By exploring the complete phase diagram of the equilibrated 
wave, we discover that for low but nonzero values of the disorder strength, a large-scale structure 
-a condensate- appears in the equilibrium distribution. 


Because of their high-degree of control and tunabil- 
ity, quantum gases have become versatile model systems 
to study effects originating from many different fields of 
physics such as condensed matter, quantum information, 
qiiantum hydrodynamics and even high-energy physics 
[J-Q . The plwsics of disordered systems does not escape 
the trend [3,3 with the observation of coherent backscat- 
tering and Anderson localization with non-interacting 
matter waves UM- When interactions are addition¬ 
ally present, disordered gases offer even richer phenom¬ 
ena. For instance, the phase diagram of interacting dis¬ 
ordered Bose gases (“dirty boson” problem) at zero and 
finite temperatures has recently stirred considerable the¬ 
oretical and experimental interest [Tsl - f^ . 

Another important, yet poorly understood problem, is 
the long-time limit of the out-of equilibrium dynamics of 
quantum gases released in a disordered potential. For 
spatially-narrow initial states, the atomic cloud spreads 
out and, in the absence of interactions and under suit¬ 
able conditions, eventually freezes due to Anderson lo¬ 
calization [ 2 ^, [^. However, if the matter wave is - 
even weakly- interacting, theoretical and experimental 
evidence suggest that the cloud keeps expandin g in defi- 
nitely so that no stationary state is ever reached (2514^ . 
When atoms are prepared in a plane-wave state at some 
finite velocity, the interesting dynamics takes place in mo¬ 
mentum space [l3. [28l - [^ and the nature of the system at 
very long times, resulting from the complicated interplay 
between atomic collisions and scattering off disorder, is 
largely unknown. This is the question we address in the 
present Letter. 

Intuitively, the dynamics of weakly-interacting matter 
waves in momentum space is expected to yield a thermal¬ 
ization process due to atomic collisions. In the homoge¬ 
neous case, this process has been extensively studied for 
nonlinear waves obeying the nonlinear Schrodinger (NLS) 
equation, including dilute Bose gases in the mean-field 
regime, within the framework of weak turbulence the- 
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ory [ 3 ^ : the system generically equilibrates to a thermal 
Rayleigh-Jeans distribution maximizing entropy [H, [s^. 
For inhomogeneous systems like disordered systems, this 
thermalization process needs further analysis since the 
density of states (DoS) is dramatically altered, in particu¬ 
lar at low energies. We have found that, if the equilibrium 
momentum distribution achieved by a weakly-interacting 
matter wave evolving in a two-dimensional (2D) speckle 
potential is still a thermal Rayleigh-Jeans distribution, it 
also exhibits, for specific values of the disorder strength 
and of the initial velocity, a large-scale coherent struc¬ 
ture. This “condensate” coexists with the background of 
thermalized atoms and disappears at vanishing disorder. 
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FIG. 1. Cut at fcy = 0 of the equilibrium momentum distribu¬ 
tion Uk/N, obtained numerically for Vo = Q.75eg, Eq = 1.5£<; 
and number density n = 2/^^. Dashed curve; <7 = 0, see Eq. 
©• Solid curve; 5 A 0, see Eq. ®. Insets; same distribu¬ 
tions shown in the plane (kx,ky) for g = Q (up) and g yf 0 
(down). 

Quench scenario - We consider a gas of N weakly- 
interacting bosons, initially prepared in the plane-wave 
state Ifco), and suddenly released at t = 0 in a 2D random 
potential VYr) that we chose to be a blue-detuned speckle 
potential [7|. Such a random potential has equal mean 
and root-mean square amplitudes Vq and we use its corre¬ 
lation length C and its correlation energy 
as units of length and energy. After the quench, the 
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cloud expands according to the NLS equation, ihdt'i' = 
—/i^/(2m)V^'I' + V{r)'i> + (/I'f'pvf' with j'l'p normalized 
to the total number of atoms N. Finally, after a time 
t, the speckle potential is switched off and the disord er- 
averaged momentum distribution nk(t) = of 

the cloud is recorded. We are interested in the long-time 
asymptotics —>• oo) = In the absence of interac¬ 

tions ((7 = 0), this limit is typically established after the 
momentum distribution is isotropized by scattering, i. e. 
after a few transport times r, and reads [28l - l3lj| : 

f J / n\ n\ 


Here A^{k) is the disorder-averaged spectral function of 
the speckle potential and ^(e) = / k/A^{k) is the 
disorder-averaged DoS per unit volume, see Supplemen¬ 
tary Material. Since / deAe{k) = 1, we see that Ae{ko) 
is the probability density that an atom with momentum 
fco has energy e after the quench. Scattering from the 
speckle being elastic, this energy distribution remains 
constant during the evolution and, finally, when the po¬ 
tential is switched off, atoms with energy e acquire a 
momentum k with probability density A^{k)/v{e). Note 
that Eq. © is in fact only valid away from the direc¬ 
tions ±ko, where coherent back and forward scattering 
interference peaks are expected (l^. [28l - [^ . When g 7 ^ 0, 
these peaks turn out to disappear over a time scale much 
shorter than the thermalization process discussed in the 
paper. We will thus not consider them in the following, 
leaving their study in the presence of interactions for fu¬ 
ture work. For an isotropic speckle, A^(k) only depends 
on |fc| = k and Uk is isotropic as well. The momentum 
distribution ([T]) is shown in Fig. [T]for Vq = 0.75e^ and 
Eq = H^Uq/ {2m) = 1.5e,^. It has a characteristic ring 
shape reflecting the energy profile of the spectral func¬ 
tion [ 2 ^. 

Thermalization ~ When g ^ 0, atoms are scattered 
both from the fluctuations of V (r) and of the nonlinear 
(random) potential gl'i'{r, t)\^■ The latter process redis¬ 
tributes energies over a collision time scale denoted by 
Tcoii. In what follows, we assume Tcoii ^ t, which is 
fulfilled for small enough g. Physically, this condition 
means that scattering events on the random potential 
occur more frequently than atomic collisions. Conse¬ 
quently, the disorder isotropizes atomic momenta before 
the nonlinearity starts to play a role. Thus, when g ^ 0 
the momentum distribution at t ^ r is still isotropic, 
but it keeps evolving: 


nk{t » r) 
N 


= J deA^{k) 


/(c<) 


( 2 ) 


which generalizes Eq. ©. The energy distribution /(e, t) 
of the matter wave is normalized according to n = N/fl = 
J de z^(e)/(e, t), where il is the volume of the system. The 
change in time of /(e, t), due to atomic collisions, is con¬ 


trolled by the kinetic equation 

= J dundu,de,W,,,,.,3.,, [{feUh, 

+ fe2feje^ - /e/es/ea “ fefe2fej , (3) 

where we have used the shorthand notation /g. = 

f{ei,t) and where k^. = y/2mei Ui/h (u^ is a unit 
vector). The collision kernel is = 47r x 

{2Trg/h)'^iz{e2)i'{e3)’z{e4)d{k^+ke^-k^^-k^J6{e+€2-e3- 
64 ), where the two delta functions stem from momentum 
and energy conservation during a collision process. Eq. 
© describes the effect of atomic collisions in a “disor¬ 
dered background” and reduces to the usual kinetic equa¬ 
tion describing short-range interactions between bosons 
[ 3 ^ in the homogeneous case where jz(e) = m/{ 2 'Kh?). 
In the absence of interactions, df{e,t)/dt = 0, whence 
f{e,t) = = O'*") = nA^{ko)/v{e) and one recovers 

Eq. ID). Furthermore, Eq. |D guarantees energy conser¬ 
vation at all times. Expressing it immediately after the 
quench, and assuming the interaction energy is negligi¬ 
ble, we find Etot = f de eiz{e)f{e, t) = N{Vo + Eq). 

The equilibrium energy distribution is obtained 
by canceling the collision kernel and thus solving for 
+ /£2/e3/£4 - fefeje^ “ hf^2!^A = 0 ^ith £ = 
£3 + €4 — € 2 - One finds the Rayleigh-Jeans distribution, 
feq _ equilibrium momentum distribu¬ 

tion thus reads: 


'Ht - f r] ^e(fc) T 

N J ^ n e— g 


{9 ^ 0 ), 


( 4 ) 


which replaces Eq. ©• The blue-detuned speckle po¬ 
tential being bounded from below by zero [39|, energies 
start from 0 in the integral. Since > 0, we must have 
T > 0 and g < min(e) = 0. The parameters T and g, can 
thus naturally be interpreted as the equilibrium temper¬ 
ature and chemical potential of the cloud after its out- 
of-equilibrium evolution. They are obtained from atom 
number and total energy conservation: 


rSmax / \ 

n = T de—and 

Jo e - M 

n{Vo + Eo)=T dee^. 

Jo e - M 


( 5 ) 

( 6 ) 


As i'{e —>• oo) = ml{2TThJ), both integrals have an 
ultraviolet (UV) divergence and an upper cutoff Emax 
has been introduced for regularization. This divergence 
stems from the mean-field description of interactions, en¬ 
capsulated by the NLS equation, which provides a clas¬ 
sical nonlinear wave model of the gas only valid at low 
energies. Following 32|, we proceed by analogy with a 
non-interacting Bose gas at equilibrium and use the cut¬ 
off Ejnax = T -\- g, which physically describes the cross¬ 
over region separating the classical-field and particle-like 
(Boltzmann) descriptions of the gas. 
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Once Eqs. (O and dG]) have been numerically solved for 
T and fx, for given values oi Eq, Vq and n, the thermalized 
momentum distribution follows from Eq. dH) . The result 
is shown in Fig. [T]for Eq = 1.5e^, Vq = 0.75e^ and n = 
2/C^. Unlike for g = 0, the distribution is now centered 
at /c = 0 and has a smooth shape, (except in the far tails, 
whose details depend on the UV cutoff). 



FIG. 2. Main panel: equilibrium phase diagram {Eo,T) ob¬ 
tained from Eqs. m and for various values of Vo at num¬ 
ber density n = 2/('^. We find a critical curve (black dots) 
separating a thermal region (blue region, < 0) from a re¬ 
gion where condensation occurs (red region, = 0). Solid 
curves: temperature T as a function of Eq in the thermal re¬ 
gion. Dashed curves: temperature T as a function of Eq in 
the region where condensation occurs. The value of Vo is in¬ 
dicated on top of each equipotential. The lowest black curve 
(Vo = 0) reproduces the scenario of weak turbulence in two 
dimensions, for which there is no condensation at finite tem¬ 
perature (Tc = 0). Inset: chemical potential /r as a function 
of Eq, for Vo = 0 and 1.2e^. When Vo 0, jx vanishes at 
Eo = Ec- 

Condensation - We now explore in more details the 
(T, iilo)-phase diagram of the system for different disor¬ 
der strengths Vq at fixed number density. Fig. [2] shows 
our results obtained by solving numerically Eqs. m and 
for several values of Vq and Eq at n = 2/ For values 
of Vq and Eq in the red region, no solution fulfilling the 
constraints T > 0 and /i < 0 can be found. The bound¬ 
ary of this red region (black dots) defines a critical curve 
Ec{Vo,n) where the chemical potential fj, vanishes (the 
inset shows the vanishing of fx for Vq = l-2e(j). To each 
value of Ec corresponds in turn a certain critical temper¬ 
ature Tc(Vo,n). Physically, this critical line signals that 
at fixed (Vb, n), there is a saturation of the population of 
excited energy levels, which can no longer accommodate 
particles if Eq is decreased below Ec (and thus when T 
becomes smaller than Tc). This phenomenon is charac¬ 
teristic of the appearance of a Bose-Einstein condensate 


(BEG) in the ground state of the speckle potential [4l|. 
For T < Tc, the gas thus consists of a thermal part co¬ 
existing with A^bec condensed atoms. In this phase, the 
total energy of the system is still given by Eq. dH]), but 
with /r = 0 as in the usual Bose condensation. Solving 
this equation then gives access to the temperature for 
Eo < Ec and is shown in Fig. [5] (dashed curves). Let us 
briefly discuss the shape of this low-temperature phase. 
First, condensation in the speckle potential exists only 
at low enough values of Eq and Vq. The reason is that 
at large Eo or Vb, too much energy is injected in the sys¬ 
tem at t = 0, and the final equilibrium temperature is 
correspondingly too large for condensation to be possi¬ 
ble. Second, it should be noted that the presence of the 
random potential is crucial for the emergence of a BEG. 
Indeed, in the limit of vanishing disorder, obtained by 
setting Vb = 0 and v{e) = m/{2'Kh'^) in Eqs. ® and 
/r vanishes only for Eq = 0 (see the inset of Fig. ®, 
and Tc —>■ 0: no BEG ever appears. In fact, this limit 
coincides with the scenario of weak turbulence of nonlin¬ 
ear waves [ 3 ^, for which it is well known that there is 
no BEG formation in 2D homogeneous systems except at 
T=0 [ 3 ^. The reason for this difference lies in the be¬ 
havior of the DoS, which vanishes at low energies in the 
presence of disorder, while it is always constant in the 
absence of disorder. 

To find the fraction of condensed atoms when T < Tc, 
we must update Eq. ([4]) by adding the BEG component: 


nk ciBEci^) 
iV “ N 


Jo ri e 


{T < Tc), 


(7) 


where nBEc(fe) is the condensate momentum distribution 
with J d^A;/(27r)^ nBEc(fc) = A^bec- Integrating over k, 
we get the fraction of condensed atoms Nbec/N = I — 
(T/n) fJ" d€iy(e)/e. Once T is known from Eq. ([5]) with 
fx = 0, one can use Eq. (0 to find Nbec/N, see Fig. [31 



FIG. 3. Condensate fraction as a function of Eq (in units of 
the critical energy Ec{Vo,n) where condensation occurs) for 
several values of Vo- The total number density is n = 

Inset: comparison between the condensate fraction at Vo = 
0.3e,^ and the theoretical prediction ([8]) valid near Eq = Ec 
(dashed line). 
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We see that the condensate fraction is always smaller 
than unity: because of the quench at t = 0 , the average 
energy Vq is suddenly injected to the gas, which eventu¬ 
ally leads to a rather high equilibrium temperature. By 
expanding the condensate fraction around T = and 
eliminating T, we obtain the following theoretical pre¬ 
diction valid for Eq close to £'c(Vb,n): 

.^BEC Ec — Eo I + Tci'{Ec)/n 

“ 1 ^ “ Ee 1 + Vo/E, + Tlv(E,)l(nE^)' 

Prediction ([51) is compared to the exact result at Vb = 
O.Se,; and n = 2/C^ in the inset of Fig. |3l 

To evaluate the condensate momentum distribution 
n-BEcik) and its typical width Afc, and further find 
when T < Tc, we proceed by analogy with the study of 
the ground state of a Bose gas at equilibrium in a confin¬ 
ing potential El- Condensation takes place in the low¬ 
est energy state of the system which is found, since the 
blue-detuned speckle potential is bounded from below, 
in the largest potential well, of typical size R ~ 1 /AA:. 
If the system size -x/fl is large enough, the condensed 
state lies in the Lifshitz tail of the density of states 
where i^(e) ~ exp(—e^/e) [d^, and its energy is typ¬ 
ically e{R) ~ hF'/(2mR^). The radius R of the con¬ 
densed state is such that il deiy(e) = 1. This gives 
Afc ~ 1/R ~ ln“^^^(-\/f2/C)- From the normaliza¬ 

tion of nBEc(fc), we have nBEc(fe = 0 )(Afc )2 ^ Abec, 
such that nBEc(fc = 0 ) ^ ln(-\/ri/C)- Thus, in the 

thermodynamic limit —>■ oo, we see that the momen¬ 
tum distribution for T < Tc consists of a narrow peak 
at fc = 0 , riBEcik) —>■ Abec^(^)i sitting on top of the 
Rayleigh-Jeans background of thermal atoms. 

Time scale of thermalization - We now address the 
question of the time scale TcoIi needed to achieve ther¬ 
mal equilibrium. To this end we consider a small per¬ 
turbation /(.(t) = ff^ -\- A(<)(5(e) and look at its expo¬ 


nential relaxation A{t) = A(0) exp (—t/rcoii) when sub¬ 
stituted in the linearized kinetic equation obtained from 
©• This procedure, detailed in the Supplementary Ma¬ 
terial, leads to Tcoh/t^ = [ti^/{gm)YE{EQ,Vo,n), where 
tq = mCf^/h. The dimensionless function E is smooth as 
long as one is not too close to the condensation threshold. 
For Eq = 1.5e,^, Vq = 0.2e^ and n = 2/^^, the gas is in the 
thermal phase, see Fig. [21 and we find F{Eo,Vo, n) ~ 0.6. 
For a quasi-2D gas of atoms with a t ypi cal s-wave 
scattering length a = 200ao (oq = 53pm) [2^ and con¬ 
fined in the transverse direction by a harmonic potential 
of frequency = lOOHz, we find Tcoh — 550t^. Tak¬ 
ing the value f = 0.2/rm, achieved in recent experiments 
111 ], we obtain TcoIi — 14 ms, a value much smaller that 
the typical duration of experiments (of the order of a few 
seconds). This makes the observation of the thermaliza¬ 
tion process described in this Letter accessible to current 
experimental setups. 

Conclusion - We have shown that the momentum 
distribution of a nonlinear matter wave released with 
finite velocity in a 2D speckle potential undergoes a 
thermalization process to a Rayleigh-Jeans distribution. 
At low enough (but nonzero) kinetic energy and disor¬ 
der strength, a condensation process also takes place at 
fc = 0. These predictions should be within reach of cold 
atoms experiments. Future studies could include a de¬ 
tailed analysis of the dynamics of thermalization and con¬ 
densation, and address the effect of residual interactions 
in the equilibrated state, which for static systems at equi¬ 
librium lead to the BKT physics [T^. Finally, it would 
be interesting to consider the case where the initial state 
is no longer a plane wave but a wave packet. This would 
allow to establish a connection between the dynamics in 
momentum space and in configuration space |25l - l27l | . 
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Appendix A: Spectral function and density of states 
in a 2D speckle potential 

In our numerical simulations, I^(r) is chosen to be 
a 2D , blue-detuned speckle potential with mean va lue 
V^(r) = Vo and two-point correlator (5y(r)5y(r') = 
[2VoJi{\r - r'|/C)/(|r - r'|/C)]^ where Jt/(r) = y(r) - 
Vo and is the correlation length. This poten¬ 
tial is numerically generated by convoluting a circu¬ 
lar Gaussian random field with a cutoff function that 
simulates the diffusive plate used in experiments [I|, 
1^. The disorder-averaged spectral function, A,.{ko) = 

J{dt/h){ko\exp{—iHt/h)\ko)exTp{iet/h), is obtained by 
propagating the plane wave state |fco) with the Hamil¬ 
tonian H = p^/(2m) -I- V, taking the Fourier transform 
of the result with respect to time and redoing the same 
calculation for many disorder configurations (as indicated 
by the over bar) [3|. The time propagation is achieved 
on a 2D grid of size L x L with periodic boundary condi¬ 
tions along X and y, by using an iterative method based 
on the expansion of the evolution operator exp{—iHt/h) 
in combinations of Chebyshev polynomials of the Hamil¬ 
tonian ia. In the simulations, a cell of surface (ttC)^ 



FIG. 1. 2D density of states i^(e) in a speckle potential, for 
Vo = O.fSe^. Inset: spectral function ^^(feo), for Vo = 0.75e,; 
and |feo| = 1.7/G 
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is discretized in typically 8-10 steps along both axis x 
and y and results are averaged over 8 disorder realiza¬ 
tions. The inset of Fig. [T] displays a typical plot of the 
disorder-averaged spectral function Ae(fco) as a function 
of the energy e, for Vo = O.TSe^ and |fco| = 1-7/C with 
e( = the correlation energy of the speckle po¬ 

tential. 

Finally, the disorder-averaged density of states per unit 
volume in the speckle potential, jz(e), is obtained by com¬ 
puting spectral functions at different values of k, and 
then summing then over k. Fig. [T] shows ^{e) obtained 
with this procedure for Vo = 0.75e^. As expected, ix{e) is 
zero for e < 0 (the blue-detuned speckle potential is al¬ 
ways positive), and iz{e) converges to the free-space result 
m/{2'Kh^) at large energies. 


Appendix B: Estimation of the thermalization time 

From the kinetic equation (3), we can estimate the 
collision time Tcou, which gives a lower bound for the 
time needed by the system to reach the Rayleigh-Jeans 
thermal distribution = T/{e — y) when g ^ 0- For 
this purpose, we consider a small perturbation from the 
equilibrium distribution at a certain energy eg: 

/,(t) = /r + A(t)5(e-eo). (Bl) 

Substituting Eq. (ED for all the in Eq. (3) and 
linearizing in A(t), we obtain, after performing integrals 
over u, U2, Ug and U4: 

- eo) = Sg'^hA{t) j de2de3de4 

+ e .2 - £3 - (^4)xike,ke2,ke3,ke4) 

x{(5(e4 - eo) [Uife + ^ 2 ) - fje 2 ] + 

<^(£3 — £ 0 ) [feiife + /ea) — /e/e 2 ] + 

^(£ — £ 0 ) [fesfei — /e2(/e3 + /e4)] + 

<^(£2-£o)[/.3/.4-/.(/.3+/.4)]}- (B2) 

The function y has been explicitly calculated in Q. For 
clarity we do not reproduce its expression here. After 
integration of both sides of Eq. (IB2p over an infinitesimal 
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FIG. 2. Plot of F{Eq, Vo, n) as a function of Eq, for Vb = 0.2e,; 
and n = 2/^^. 


energy interval centered around eoi only the third term 
in the right-hand side survives, which yields 


A{t) oc exp 


t 


(B3) 


with 


1 

T'coiKeo) 


Sg'^h 


pT 


de2d€3de4i^(€2)i'(€3)iy(e4) 


xS{eo + e 2 - €3- e4)x{keo,ke2,ke3,hi) 
x[feJe,-fe2{U+feJ], (B4) 


where we have used the same ultraviolet cutoff as in the 
main article. Eq. (IB4I) is the general expression for the 
collision time as a function of the energy cq. Since we 
only aim at giving an estimation of Tcoh, we focus on 
the particular energy eo = 0 (where the Rayleigh-Jeans 
equilibrium energy distribution is maximum). Explicit- 
ing X i and performing the integral over € 2 , Eq. (Ell 
simplifies to 


7-0011(0) 



2 

F(Eo,Vo, n), 


(B5) 


where = mC,'^/h. The inverse of the function F is given 
by: 


F~^{Eo,Vo,n) = dx dyi>{x + y)i>{x)i>{y) x 

Jo Jo 

-47Tp.f‘^Q{y/x + y/y- ^/x + y)Q{y/x + y- - y/y\) 
y/¥y{x - fj.){y - y){x + y - y) 


(B6) 


where 0 is the Heaviside function and where we have 
introduced the dimensionless density of states z>(e) = 
fi^i/(e)/m, temperature T = T/e^ and chemical poten¬ 
tial y = y/cQ. 

Fig. [2] displays F{Eo,Vo,n) as a function of Eq, for 
Vq = 0.2e(j and n = E increases with Eq at 

large enough Eq, which is consistent with the intuitive 
picture that it is harder to thermalize when more en¬ 
ergy is injected in the system. The sudden increase of 
F at small Eq signals the approach to the BEG transi¬ 
tion, which for Vq = 0.2e^ and n = takes place at 
Eq = Ec — 0.9e^ (see Fig. 2 of the main article). Note 
that the knowledge of F for E < E^ would require a more 
general kinetic equation describing exchanges between 
condensed and non-condensed particles, a task which is 
beyond the scope of the present work. From Fig. [H we 
find E{Eq, Vq, n) ~ 0.6 for Eq = l.Se^, which is the value 
used in the main article. 

For completeness we mention that F oc at low 

density, a property that we have checked numerically. 
Making use of Eq. (IB5I) , we thus find that the condition 
Icon ^ T, expressing the validity of our kinetic approach, 
reduces to lh^/(gmn(^)]^ :$> t/tq. The speckle scattering 
time T has been calculated in [j for weak disorder, where 
it was found that t/tq ~ (ec/^)^ energies. With 

this expression, and overlooking the Vb-dependence of F, 
the condition of weak interactions quite naturally reads 
gn <C Vq. 
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